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We devise an inverse statistical-mechanical methodology to find optimized interaction potentials 
that lead spontaneously to a target many-particle configuration. Target structures can possess 
varying degrees of disorder, thus extending the traditional idea of self-assembly to incorporate 
both amorphous and crystalline structures as well as quasicrystals. For illustration purposes, our 
computational technique is applied to yield an optimized isotropic (non-directional) pair potential 
that spontaneously yields the three-coordinated honeycomb lattice as the ground state structure 
in two dimensions. This target choice is motivated by its three-dimensional analog, the diamond 
lattice, which is known to possess desirable photonic bandgap properties. 

PACS numbers: 82.70.Dd, 81.16.Dn 



"Self-assembly" typically describes processes in which 
entities (atoms, molecules, aggregates of molecules, etc.) 
spontaneously arrange themselves into a larger ordered 
and functioning structure. Biology offers wonderful ex- 
amples, including the spontaneous formation of the DNA 
double helix from two complementary oligonucleotide 
chains, the formation of lipid bilayers to produce mem- 
branes, and the folding of proteins into a biologically ac- 
tive state. Molecular self-assembly is a potentially power- 
ful method to fabricate atomicallyprecise materials and 
devices. For example, Whitesides [l| has shown intricate 
two-dimensional patterns can emerge in self-assembly of 
organic molecules on an inorganic surface. Jenekhe and 
Chen have devised 'smart plastics' that assemble into 
photonic crystals. Manoharan et. al. have self-assembled 
unique, small clusters of microspheres These exam- 
ples provide glimpses into the materials science of the 
future, i.e., devising building blocks with specific inter- 
actions that can self-organize on a set of larger length 
scales. Theoretical work has mainly focused on finding 
the structure and macroscopic properties of many-body 
systems given the interactions - what we refer to as the 
"forward" problem of statistical mechanics. The forward 
problem has been extensively studied in the context of 
the freezing transition both analytically 14 and numer- 
ically and more recently by Kamien '6!|, who uses 
geometric arguments to obtain crystal entropy. 

The purpose of this Letter is to introduce an inverse 
statistical-mechanical methodology to find optimized in- 
teraction potentials that lead spontaneously to a tar- 
get many-particle configuration. The so-called 'reverse' 
Monte Carlo method 7, 8J has been used to obtain in- 
teractions in liquids given the pair correlation function, 
which only has partial configurational information. Our 
inverse methodology distinguishes itself in that we apply 
it to self-assembly of a given N-particle configuration, 
which may be be crystalline, quasicrystalline, or amor- 



phous. We envision target structures possessing varying 
degrees of disorder, which enables us to extend the tra- 
ditional idea of self-assembly. 

The idea of tailoring potentials to generate targeted 
structures is motivated by the rich array of fundamen- 
tal issues and questions offered by this fascinating in- 
verse statistical-mechanical problem as well as our re- 
cent ability to identify the structures that have optimal 
or desirable bulk properties. The latter includes novel 
crystal structures for photonic band-gap applications j|| , 
materials with negative or vanishing thermal expansion 
coefficients ^|, materials with negative Poisson ratios 
jllj . materials with optimal transport and mechanical 
properties [l^ , and mesoporous solids for applications in 
catalysis, separations, sensors and electronics [l^ . 

Colloids are the ideal system to test our optimized po- 
tentials, since both repulsive and attractive interactions 
can be manipulated (e.g., via depletion forces, dipole- 
dipole interactions, electrostatic interactions, etc. Q) 
and therefore offer a panoply of possible potentials that 
far extends the range offered by molecular systems. 

Because there is a vast class of many-body potentials, 
we will focus on isotropic pairwise additive interactions 
for simplicity here. There are many open questions even 
for this simple class of potentials. For example, it is 
not known what are the limitations of isotropic pair- 
wise additivity for producing target structures. We know 
that such interactions cannot produce thermodynami- 
cally stable chiral structures with a specified handedness; 
equal amounts of left-handed and right-handed struc- 
tures would result. When is anisotropy in the potential 
required? An answer based on intuition from molecu- 
lar systems would fail here. For instance, the diamond 
lattice is thought to require directional interactions be- 
cause such structures found in Nature result from cova- 
lent bonding. In fact, it is not known whether a diamond 
lattice could be created from an isotropic pair potential. 
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This structure has a special status in photonics research 
because a diamond lattice of dielectric spheres exhibits a 
photonic band gap across the Brillouin zone 0. 

The two-dimensional analog of this open three- 
dimensional crystal is the three-coordinated honeycomb 
lattice. Accordingly, our general optimization procedure 
(described below) will be illustrated by applying it to 
produce an optimized circularly symmetric pair potential 
V{r) that spontaneously yields the honeycomb lattice as 
the ground state (zero-temperature) structure in a pos- 
itive density range. In contrast to previous approaches 
that have claimed to produce open lattice structures, our 
procedure incorporates the phonon spectra, which is a 
crucial ingredient. Because the honeycomb is an open 
lattice that is a subset of the triangular lattice, it is inher- 
ently challenging to assemble using isotropic potentials. 
Indeed, such a potential has never been found before. 

The potential energy for a system of N classically in- 
teracting particles at positions r'^ = ri , r2 , . . . , tn in the 
absence of an external field is given by 

$(r^)=^y2(r„r,)+ ^ y3(r., r„ r,.) + ... (1) 

i<j i<j<k 

where the y„'s are n-body potentials. In this study, we 
consider only isotropic pair potentials and therefore 



$(rN) 



i<j 



(2) 



A central feature of our inverse approach is a com- 
putational algorithm that searches for and optimizes a 
functional form for V{r) that leads to self-assembly of a 
given target structure. To find an optimized V{r) for a 
given target structure, we make an initial guess for that 
function. We require that this initial potential have real 
frequencies for each of its normal modes (for a lattice, 
this means real phonon frequencies for all wavevectors in 
the Brillouin zone). Thus, the structure is mechanically 
stable at zero temperature. We then parameterize the po- 
tential, establishing a family of functions V{r; {oq. ..«„}) 
of which our initial guess is a member. The parameteri- 
zation must be chosen so that an overall rescaling is not 
possible. For each a^, we choose a range of values that 
it can take, namely [a™*",a™°^]. We then optimize this 
family of functions for self-assembly. Specifically, the pro- 
gram runs a molecular dynamics simulation (MD) at vol- 
ume (or area) per particle a, initially in the target struc- 
ture configuration and interacting via the initial guess 
potential. The initial root mean square speed in the MD 
corresponds to about 90-95% of the melting temperature, 
where the system is in a non-harmonic regime and liquid 
nucleation and/or structural phase transition is begin- 
ning to set in. The root mean square deviation from the 
target structure, defined as 
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is computed and averaged over a number of simulations. 
Here is the position of the i*'' particle after an ap- 
propriate amount of simulation time, rf^"^ is its initial 
position, and N is the number of particles. The quan- 
tity L is thus minimized in parameter space by simulated 
annealing. The program outputs the set of parameter 
values that corresponds to the minimum value of L, pre- 
sumeably giving the potential that best suppresses liquid 
nucleation and/or a possible structural phase transition. 
This is the idea of the algorithm: we postulate that if the 
potential is modified in such a way that deformations of 
the target structure are suppressed near (but below) the 
phase coexistence region, the structure will self-assemble 
from a random configuration in a MC simulation. In the 
case that the melting temperature rises significantly over 
the course of the optimization, the initial temperature 
can be increased such that the system again approaches 
coexistence. A more detailed description of the optimiza- 
tion procedure will be given in |l5| . 

We now apply our methodology to obtain an optimized 
circularly symmetric interparticle potential, V(r), that 
will spontaneously favor the self-assembly of randomly 
placed particles in two dimensions into the honeycomb 
lattice upon simulated annealing from high to zero tem- 
perature. A claim was made that the honeycomb lattice 
was favored by a hard-core plus linear-ramp potential for 
certain parameter values [l6j, but we have now shown 
conclusively that this is not possible because the phonon 
spectra of such V{r) have imaginary frequencies. Thus, it 
would be a significant accomplishment for our procedure 
if it assembled the honeycomb lattice. 

For lattice self-assembly, we make another restriction 
on the initial function: that the target should be energet- 
ically favored among the four principal 2D lattices, (the 
triangular, square, honeycomb and Kagome) over a sig- 
nificant range of area per particle a. This is the second 
of our two necessary conditions for lattice self assembly, 
the first being that all phonon frequencies are real. 

Both the triangular and honeycomb lattices have their 
first near neighbors (in relative distances) at 1, -x/S and 2. 
The coordination numbers for those neighbors are 6, 6, 
6 and 3, 6, 3 for the triangular and honeycomb lattices, 
respectively. Thus, we choose a potential V{r) such that 
V{1) is not negative, or else there will be a tendency to 
accumulate as many neighbors as a lattice will permit, 
and hence it will fall into the triangular. Yet if there is 
no local minimum at the nearest neighbor distance, the 
mechanical stability (phonons) will almost certainly be 
removed (as we have found after trial and error). This 
motivates the choice of the family of functions 
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This is not the function that is finally given to the op- 
timization program; here, we are still approximating the 
potential form. This is just a LJ-type interaction (12/10 
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FIG. 1: V(r), potential that favors honeycomb lattice self- 
assembly. 
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FIG. 2: Phonon spectrum (frequency squared) for honey- 
comb with potential V{r) at a = 1.45. 



rather than 12/6) added to a soft repulsive exponential 
tail that is normalized to 2TrA in 2D regardless of the 
value of A. We choose the (12/10) LJ because it is more 
sharply distinguished from the exponential tail than is 
the (12/6). It also has its minimum at unity and has 




FIG. 3: 500-particle annealed MC results, for potential with 
parameters displaced from initial guess, a = 1.45. 
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FIG. 4: 500-particle MC results, annealed from knT = 0.22 
to fcsT = at a = 1.45. 

unit depth. For certain choices of A and A, this potential 
has one local minimum but with a positive value. To find 
the best values for A and A to stabilize the honeycomb, 
we choose an a at which the honeycomb nearest neigh- 
bor is at unit distance and we maximize the interstitial 
hopping energy (both local and distant) of just the ex- 
ponential term over A at a given A. We can then set a 
lower bound on the value of A by constraining the po- 
tential to have a positive hopping energy (hopping into 
honeycomb interstitial sites is unfavorable). A qualita- 
tive upper bound is set on A just by the need for mechan- 
ical stability: for large enough values of A, this potential 
becomes purely repulsive and the honeycomb is mechan- 
ically unstable with it. So A is chosen to be its lower 
bound. This procedure gave the values A — 3.0 and 
A = 2.677. The resulting potential gave favorable lattice 
sums but still imaginary phonon frequencies, so in order 
to 'brace' the lattice, an extra attractive Gaussian was 
added, and so we use the paramterization 

V{r;ao,ai,a2,a3) = A_ ^+a^e-'^-'^-0.4e-40('— 3)^ 

(5) 

For parameter values oq = 6.0, ai — 21.5, 02 = 2.677, 
and 03 = 1.829, the guess potential meets our two neces- 
sary conditions. A 500-particle annealed MC simulation 
using this potential produced a lattice reminiscent of the 
honeycomb, but with a significant number of defects in 
the ground state. 

In order to demonstrate the effectiveness of the opti- 
mization program, we somewhat arbitrarily displaced the 
parameters from the initial guess function, setting them 
to be ao = 6.5, ai = 18.5, 02 = 2.45, ag = 1.83. The 
resulting 500-particle annealed configuration is shown in 
Fig. 131 We then started the optimization program with 
these values for the parameters, the output of which was 
the following potential: 

Vf{r) = A-^ + 17.9e-2■49-_o.4e-4"(--l■823)^ (e) 
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This function is plotted in Fig. ^ The phonon spec- 
trum is given in Fig. |21 The lattice sums (not shown 
there) demonstrate that while the honeycomb is stable 
over a wide range of a, the triangular eventually dips 
below it. The global minimum for the triangular occurs 
at the area at which the first nearest neighbor lies at the 
bottom of the Gaussian. As long as there is an attractive 
gaussian in V{r), there is no avoiding this effect. How- 
ever, it has little relevance for the lattice self-assembly, 
since it is at relatively high a. 

Monte Carlo simulated annealing was carried out on 
a 500-particle system (at a = 1.45) interacting via the 
potential shown in Fig. (equation EJ . The resulting 
configuration is depicted Fig. ^ Except for a few de- 
fects, the honeycomb has indeed self-assembled. The de- 
fects actually seem to be "missing particles," rather than 
dislocations or pockets of disordering. One might expect 
that increasing the number of particles to fill the defects 
(or adjusting a accordingly) would eliminate the defects, 
but this is not the case. The defects are likely due to the 
slow dynamics of the MC, i.e., they were "frozen in" dur- 
ing annealing. According to the lattice sums, the perfect 
honeycomb lattice is lower in energy than the one pro- 
duced in the MC simulation with defects. 

We have found that as long as the salient features of 
the honeycomb potential are kept (two local minima at 
distance ratio -s/S, the first being positive and the second 
negative), self-assembly is unaffected by perturbations in 
the potential, i.e. the potential is robust. This is essential 
if this system is to be tested experimentally. 

In summary, using an inverse statistical-mechanical ap- 
proach, we have found an optimized isotropic pair po- 
tential that results in the self-assembly of the targeted 
honeycomb lattice. In many nanoscopic systems, exper- 
imentalists have increasingly greater control over inter- 
component interactions, and hence, optimal design of 
nanostructures by self-assembly ultimately is always an 
inverse problem. Our results give some hope to the pos- 
sibility of self-assembly of the more challenging diamond 
lattice with isotropic pair potentials. As we indicated ear- 
lier, this would have important implications for photonics 
devices. We are currently searching for a parameteriza- 
tion that upon optimization of an isotropic pair potential 
will yield the diamond crystal as its ground state. Note 
that our methodology also offers the opportunity to ex- 
amine defects in these ideal structures and to minimize 
their occurrence by maximizing their energy costs. 

There are many fascinating research avenues that we 
can explore using our inverse approach. Our results beg 
the question of whether more exotic structures can be 
assembled using only isotropic pair potentials. For ex- 
ample, one might try to assemble a Buckyball in a NVT 
annealing simulation with 60 particles with a potential 



that has sharp minima at the first several neighbor posi- 
tions. Although we do not know whether this is possible, 
we have already found a potential that produces small 
clusters of particles (e.g., simplices), as well as one that 
produces long chains, or "colloidal nanowires." Clearly, 
there must be limits on the types of structures that can 
be assembled using only isotropic potentials. We know, 
for example, that chiral structures with a target chiral- 
ity cannot be formed, but beyond this specific case we 
know very little about the limitations of isotropic pair 
potentials: a fundamentally important problem. 

The optimization algorithm proposed in this Letter is 
only one approach to the inverse problem, and we ex- 
pect that others will be needed to search for interac- 
tions (isotropic or not, additive or not) that stabilize non- 
periodic systems. Apart from any particular algorithm, 
however, a central point of this Letter is to propose the 
use of powerful inverse statistical-mechanical techniques 
to exquisitely control self-assembly from the nanoscopic 
to microscopic scales. 

This work was supported by the Office of Basic Energy 
Sciences, DOE, under Grant No. DE-FG02-04ER46108. 
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